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We propose a new method to obtain kinetic properties of infrequent events from molecular dy¬ 
namics simulation. The procedure employs a recently introduced variational approach [Valsson and 
Parrinello, Phys. Rev. Lett. 113, 090601 (2014)] to construct a bias potential as a function of several 
collective variables that is designed to flood only the associated free energy surface up to a prede¬ 
fined level. The resulting bias potential effectively accelerates transitions between metastable free 
energy minima while ensuring bias-free transition states, thus allowing accurate kinetic rates to be 
obtained. We test the method on a few illustrative systems for which we obtain an order of mag¬ 
nitude improvement in efficiency relative to previous approaches, and several orders of magnitude 
relative to unbiased molecular dynamics. We expect an even larger improvement in more complex 
systems. This and the ability of the variational approach to deal efficiently with a large number of 
collective variables will greatly enhance the scope of these calculations. This work is a vindication 
of the potential that the variational principle has if applied in innovative ways. 

PACS numbers: 05.10.-a, 02.70.Ns, 05.70.Ln, 87.15.H- 


Molecular dynamics (MD) simulation can provide di¬ 
rect physical insight into the time evolution of molecu¬ 
lar systems, and has become an important tool in many 
branches of science. However, the energy landscape of 
many interesting systems is characterized by numerous 
metastable states separated by large kinetic barriers [T]. 
A longstanding issue in the field of computational physics 
has been understanding phenomena which occur on such 
landscapes and involve timescales much longer than those 
accessible by traditional MD. This is arguably the most 
serious limitation of this powerful technique. While sev¬ 
eral methods have been developed to enhance the sam¬ 
pling of rare events EHIS] and obtain a static picture 
of the underlying landscape, calculating quantitative ki¬ 
netic and dynamic properties remains a challenge. 

Within the framework of transition-state theory 
(TST), one can envision adding a bias potential to en¬ 
hance barrier crossing events out of local energy minima. 
The effect of the bias on the transition rates can then 
be rescaled by considering the theory of activated pro¬ 
cesses na. This is the idea behind many accelerated dy¬ 
namics methods, such as conformational flooding nanE], 
hyperdynamics [muH], accelerated MD m and several 
other methods [2QH28] . With a well-designed bias po¬ 
tential, one can obtain substantial speed-ups relative to 
unbiased MD; however, care must be taken to ensure 
that the bias vanishes at the transition state and that 
the system remains within local minima long enough to 
accumulate local averages. In other words, transitions 
in the biased ensemble should be representative of those 
that would be eventually sampled in the unbiased case in 
the long time limit. 

Thus, the key problem in these approaches has been 


designing a bias potential that leaves the transition states 
untouched, remains simple to evaluate, yet gives acceler¬ 
ation relative to unbiased MD that scales well with the 
number of degrees of freedom. A poorly designed bias 
potential can either lead to vanishing boosts relative to 
MD or give inaccurate timescales inisn]. This is because 
the true high-dimensional potential energy surface has an 
enormously large number of low-lying saddle points, and 
the likelihood that a bias potential disturbs some of these 
saddle points in a non-trivial way becomes significant as 
the dimensionality increases naEO]. An appealing alter¬ 
nate approach is to shift the attention from dynamics on 
a potential energy surface to a free energy surface (FES), 
which is a coarse-grained representation of the system 
in terms of a few collective variables (CVs). Provided 
that the CVs distinguish between the various metastable 
states or basins, one then enhances fluctuations in these 
CVs with an appropriate bias to escape local minima. 

Along these lines, Tiwary and Parrinello have recently 
proposed a method to recover the correct dynamics from 
biased metadynamics simulations [22l[3T]. Using only a 
few relevant CVs they are able to recover correct tran¬ 
sition rates by limiting the deposition frequency of the 
time-dependent bias. Their key idea is to bias the system 
perturbatively such that the frequency of perturbation is 
between the fast intra-basin and the slow inter-basin re¬ 
laxation frequencies, thereby creating bias-free transition 
states. However, the need to reduce the bias deposition 
frequency is a severe computational handicap and forces 
one to rather lengthy calculations. In spite of this, the 
method has been used to calculate the unbinding kinet¬ 
ics of ligand/protein interactions and other applications, 
demonstrating the power of using dimensionality reduc- 
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tion j32H35] . 

In this Letter we propose a new method based on the 
efficient construction of a static bias which acts on a few 
relevant CVs and only acts below a threshold free energy, 
Fc, above which the FES is unaffected. The problem of 
estimating the free energy and constructing the bias is ac¬ 
complished within the framework of a novel variational 
approach of Valsson and Parrinello [36]. Using this vari¬ 
ational approach and including a constraint that ensures 
the bias potential will not act on the transition state, we 
obtain a static bias which accelerates transitions between 
metastable states, but makes no assumptions about the 
transition pathway or reaction coordinate. Importantly, 
no a priori assumption is made of the shape of the FES, 
which may be of arbitrary complexity. Nota bene: Flood¬ 
ing the FES, which is defined in terms of well-chosen CVs, 
is in no way the same as filling or lifting up all states be¬ 
low a certain potential energy as has been done in the 
past [laisT]. 

In the following, we assume that there is a set of 
CVs, which map a given atomistic configuration, R, 
onto a finite set of bounded coarse-grained variables, 
s = {s^(R)}. The FES is then defined up to an irrel¬ 
evant arbitrary constant as, 

F(s) = -1 log y dRS{s - s(R))e-^^(*^), (1) 

where U(R) is the interaction potential and (3 = 1/k^T. 
We seek to introduce a bias potential, U(s) that acts 
on the CVs such that the bias potential will enhance 
fluctuations out of local minima of the FES, while not 
biasing the transition state ensemble. The problem of 
finding the correct bias potential that will locally “flood” 
the FES is solved by minimizing the following functional 
of Valsson and Parrinello [36] , 

= -p log ^ /rfse-ms) + J (2) 

in which p(s) is a chosen normalized probability distri¬ 
bution. The functional, U[U] is closely related to the 
relative entropy [38] and to the Kullback-Leibler diver¬ 
gence [39]. It can be shown that U[U] is a convex func¬ 
tional and that it has a global minimum that satisfies the 
following relation, valid up to a constant, [36] . 

y(s) = -F(s)-tlogp(s), (3) 

when p{s) ^ 0 and U(s) = oc otherwise. At the min¬ 
imum the CVs will be sampled according to the target 
distribution, p(s), so we can tailor the sampling through 
the choice of p(s). This is a great advantage of the vari¬ 
ational formalism. 

We now choose a variational form for U(s) that to¬ 
gether with an appropriate choice of p{s) automatically 
ensures that the bias is zero when F(s) is greater than a 


preassigned value, Fc, relative to the minimum of F{s). 
To this effect we write V (s) in the form 

U(s) = i;(s)S(-i;(s) - F^) (4) 

where v{s) is the function that is allowed to vary in order 
to minimize 0[U]. S(x) is a sigmoidal switching func¬ 
tion such that S(x) = 1 for x ^ —oc and S(x) = 0 for 
X ^ oo. We also assume that S(x) is continuous and 
differentiable while still approximating a step function 
behavior. Within the local minimum, v{s) ~ — F(s). 
Such a functional form ensures that the bias is always 
less than Fc. A schematic depiction of the definition of 
the bias and switching function is presented in Fig. 



FIG. 1. A hypothetical one dimensional free energy surface, 
F(s) is shown in black along with the switching function, 
S(—i’(s) — Fc) (dashed blue line) for a chosen cutoff, Fc. The 
resulting bias, V(s) is depicted by the red solid line, and fills 
F(s) up to the value of Fc (green solid line). 


We note that the proof of the convexity of Q[V] in 
Ref. [36] does not make any assumptions regarding the 
exact form of the bias potential, so this property of U[U] 
is fully preserved in the presence of the switching func¬ 
tion. Of course, since we have limited the variational 
flexibility of U(s), the reconstruction of the FES from 
Eq. [^ will be only approximate. However, this is of no 
relevance here since our only purpose is to construct a 
V (s) that facilitates exit from the metastable state and 
does not affect the transition region. 

We also would like to choose the target distribution 
such that for free energies below the cutoff, Fc, the bias 
compensates almost exactly the underlying free energy 
and is zero otherwise. This can be achieved if we have 
a good estimate of the free energy surface F*(s) = F(s) 
and write p{s) as: 

^ /ds'S{F*{s') - fc)* 

This choice selects only those portions of s-space in which 
F*(s) < Fc. Of course F*(s) is not a priori known, 
but can be estimated in a self-consistent procedure in 
which we first guess F*(s), then minimize U[U] using the 
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corresponding p(s) to obtain a new V{s), which in turn 
gives a new estimate of F*(s) from F*(s) = —v{s). The 
process is iterated until a converged result is reached, in a 
way very similar to that described in Ref. [40] in a related 
but different context. One way of making this variational 
principle practical is to expand v{s) in a basis function 
set, fk{s), 

v{s) = '^ak ■ fk{s), (6) 

k 

and use the expansion coefficients, {(T/c}, as variational 
parameters isniiin]. The set of basis functions may be 
plane waves, Chebyshev polynomials, or any other suit¬ 
able basis set, depending on the nature of the CVs and 
the problem at hand. Of course any other variational 
form for v(s) can equally well be chosen. For the opti¬ 
mization of 0[R] with respect to the set of variational 
parameters {ak} we employ the stochastic optimization 
method introduced in Ref. im, as done previously in 
Refs. |36l [40]. In this work, we have employed a Fermi- 
type switching function for S(x), 

where A is a parameter with units of inverse energy that 
determines how quickly the switching function goes to 
zero. The stochastic optimization procedure is presented 
in more detail in the SM. 

After the minimization procedure, the bias, ^(s) is 
implemented as a fixed bias which speeds up the molecu¬ 
lar dynamics time by facilitating escape over free energy 
barriers. The physical time is recovered using the rela¬ 
tionship [mug, 

fT'tot '^tot 

t* At* = 

i i 

= ^MD ’ (^) 

where A^md is the MD integration time step, ^md is the 
time as measured in a biased molecular dynamics run, 
and t* is the “real” time in an unbiased simulation cor¬ 
responding to tMD- The quantity is the 

acceleration factor which is a measure of how much the 
time is boosted, and the subscript V denotes sampling 
under the biased potential. 

We now proceed with a few illustrative examples. 
First, we consider the classical case of the Cjeq Cjax 
conformational change of alanine dipeptide in vacuum, 
which can be distinguished by the two backbone dihe¬ 
dral angles, (0, The apparent barrier height for the 
forward conversion is approximately 34 kJ/mol. We bias 
both angles 0 and -0, and use a plane wave expansion, 
'r’(s) = The procedure has been implemented 

in a private development version of the PLUMED 2 plug¬ 
in m, which will be made available in the near future. 


Details of the computational methods and optimization 
procedure along with the resulting bias potentials are 
presented in the SM. To test the efficiency of the method, 
we perform a short optimization step of the bias for sev¬ 
eral different Fc values, which gives us a series of bias 
potentials of increasing strength. Subsequently, we used 
each potential as a fixed bias for 60 independent trajecto¬ 
ries all started from the same initial Cjeq configuration, 
from which we measured the first passage time (fpt) into 
the C7ax state. As a comparison we also perform well- 
tempered metadynamics (WTMetaD) with similar pa¬ 
rameters to those used in Ref. [22| (see SM). 

Fig. I^a) shows the cumulative distribution of the un¬ 
rescaled first passage times for increasing Fc values as 
compared to WTMetaD at T = 250 K. For Fc greater 
than 24 kJ/mol, the free energy flooding method drives 
barrier crossing faster than WTMetaD. This is due to 
the inferior efficiency of WTMetaD to fill the basin rela¬ 
tive to the variational approach. This lack of efficiency in 
WTMetaD is made worse by the need of depositing the 
Gaussians at a low pace. The optimization of the bias us¬ 
ing the variational approach only needs to be performed 
once for each Fc and converges rapidly (within 1 ns) to an 
effective bias. As shown in the SM, the MD time needed 
for the optimization step is negligible as compared to the 
total MD time needed for the subsequent production runs 
and thus can be ignored when considering the efficiency 
of our approach. 

Fig. Ifb) shows the fpt distributions after rescaling 
the time according to Eq. All of the transition times 
collapse onto a single distribution, which corresponds to 
the unbiased trajectory distribution. This indicates that 
even for values of Fc which are just a few /cbT below the 
apparent free energy barrier, the method recovers the 
correct distribution of passage times. Eig.[^c) shows the 
mean first passage times as a function of Fc along with 
the associated acceleration factor. While the accelera¬ 
tion factor increases by several orders of magnitude with 
increasing Fc, due to the exponential dependence of the 
rate on the barrier height, the rescaled mean first pas¬ 
sage times are constant within the error bars and are in 
quantitative agreement with the first passage time from 
unbiased MD of 23 ± 3.8 /is. We expect that a test sim¬ 
ilar to the one illustrated in Eig. |^c) will be extremely 
useful in ascertaining the reliability of the accelerated 
timescales in complex systems. Sensitivity of the rescaled 
escape rates to the extent of flooding is typically a clear- 
cut indication that the chosen collective variable is not 
equilibrating fast enough, or that the transition states 
have been corrupted m- 

As a second example we construct a model of ben- 
zophenone in water. Since the true system is adequately 
sampled by unbiased molecular dynamics we have stiff¬ 
ened the torsions around the ketone group. This stiff 
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FIG. 2. (a) Cumulative distribution of the unrescaled first 
passage times for the Creq Crax transition in alanine dipep¬ 
tide in vacuum at T = 250 K using different Fc values (From 
left to right: Fc = 30, 28, 26, 24, and 22 kJ/mol) as compared 
to WTMetaD (black), (b) Cumulative distribution of first 
passage times after time rescaling. As a guide to the eye, an 
exponential fit to the data with r = 28.2 /is is shown (black 
line), (c) Mean first passage times plotted vs. increasing Fc 
(left axis). The shaded gray region depicts the standard error 
determined from WTMetaD. Also shown is the corresponding 
acceleration factor (blue) on a logarithmic scale (right axis). 


model allows us to treat the back and forth out-of-plane 
bending motion of the two phenyl rings around the ketone 
in the presence of molecular solvent, shown in Fig. |^a), 
as a rare event. In Fig. [^b), the FES of this motion ob¬ 
tained using the variational method with dihedral angles 
01 and 02 is shown for the transition between the two 
conformations. The FES shown in Eigure |^b) is only 


a portion of the total EES in the space of the dihedral 
angles. Since we are only interested in a limited region of 
the EES, we use products of Chebyshev polynomials de¬ 
fined in the range [—1 : 1] as basis functions, which allows 
us to construct a bias confined to the region of conforma¬ 
tional space for which we are interested. Computational 
details are presented in the SM. The apparent barrier in 
Fig. ib) is 22 kJ/mol. Using a fixed value of Fc = 12 
kJ/mol we obtain a bias which partially floods both wells 
and gives a boost in the kinetics, corresponding to a new 
barrier of ^ d/cpT. In order to ensure a stable optimiza¬ 
tion, we use 4 multiple walkers during the optimization 
step. The rates of the forward and reverse transition af¬ 
ter rescaling the trajectory using Eq. are compared in 
Fig-i c) to those from unbiased simulations performed 
at different temperatures. Within the error bars, the for¬ 
ward and reverse rates are equal as expected. 


In conclusion, in this Letter we have presented an 
efficient method to obtain kinetic properties of rare 
events using a variational procedure to construct a static 
bias in free energy space. The method is general and 
easy to implement, making it suitable for a variety of 
applications. We have demonstrated the procedure on 
two model systems and have shown that at least one 
order of magnitude of efficiency can be gained relative 
to WTMetaD [22], and several orders of magnitude 
relative to unbiased MD. More impressive gains are 
to be expected in more complex systems. With this 
new method, many other computational advantages are 
possible. Since we are able to focus separately on one 
particular metastable state, we can limit our attention 
to this particular state and use collective variables 
appropriate for this state. This is a far less demanding 
requirement than finding collective variables able to 
reconstruct the whole configuration space. Eurthermore, 
we can increase the number of collective variables 
since the space to be filled is small. Approximate but 
highly dimensional forms of bias can also be used. This 
considerably alleviates the problem, sometimes vexing, 
of finding appropriate CVs. All of these considerations 
make the method preferable not only to WTMetaD, but 
to a more straightforward approach in which the global 
free energy is calculated first and only later the states 
below Fc are selected. This work is also a clear example 
of how the variational property of U[U] can be employed 
in a novel manner. 
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FIG. 3. (a) Representative conformers and the definition of 
the CVs for the stiffened benzophenone model, (b) Reference 
free energy surface (kJ/mol) for the region of configuration 
space of interest, (c) Arrhenius plot of the rate of the intercon¬ 
version from conformer A ^ B (upper) and B ^ A (lower). 
Rates obtained from biased trajectories using Fc = 12 kJ/mol 
are show with red circles and unbiased MD are shown in black 
squares. The line is a linear fit to the biased (red circle) data. 
The slope gives an estimate of 21.0 ± 0.5 kJ/mol (upper) and 
22.8 ± 0.14 kJ/mol (lower) for the barrier height which are 
both in agreement with the apparent barrier of ~ 22 kJ/mol 
obtained from enhanced sampling. 


* parrinello@phys.chem.ethz.ch 

[1] D. Wales, Energy landscapes: applications to clusters, 
biomolecules and glasses (Cambridge University Press, 
Cambridge, U.K., 2003). 

[2] G. Torrie and J. Valleau, J. Comput. Phys. 23, 187 
(1977), 

[3] T. Huber, A. E. Torda, and W. F. Gunsteren, J 
Computer-Aided Mol Des 8, 695 (1994) 

[4] B. Roux, Comput. Phys. Commun. 91, 275 (1995). 

[5] F. Wang and D. P. Landau, Phys. Rev. Lett. 86, 2050 


( 2001 ). 

[6] E. Darve and A. Pohorille, J. Chem. Phys. 115, 9169 

( 2001 ) 

[7] A. Laio and M. Parrinello, Proc. Natl. Acad. Sci. U.S.A. 
99, 12562 (2002). 

[8] U. H. E. Hansmann and L. T. Wille, Phys. Rev. Lett. 
88, 068105 (2002) 

[91 L. Maragliano and E. Vanden-Eiinden, Chem. Phys. Lett. 
426, 168 (2006) 

[10] P. Maragakis, A. van der Vaart, and M. Karpins, J. 
Phys. Chem. B 113, 4664 (2009) 

[11] J.-H. Prinz, H. Wu, M. Sarich, B. Keller, M. Senne, 
M. Held, J. D. Chodera, C. Schiitte, and F. Noe, J. 
Chem. Phys. 134, 174105 (2011). 

[12] K. J. Kohlhoff, D. Shukla, M. Lawrenz, G. R. Bowman, 
D. E. Konerding, D. Belov, R. B. Altman, and V. S. 
Pande, Nat. Chem. 6, 15 (2014). 

[13] D. Shukla, Y. Meng, B. Roux, and V. S. Pande, Nat. 
Commun. 5 (2014). 

[14] B. J. Berne, M. Borkovec, and J. E. Straub, J. Phys. 
Chem. 92, 3711 (1988). 

[15] H. Grubmiiller, Phys. Rev. E 52, 2893 (1995). 

[16] E. M. Muller, A. de Meijere, and H. Grubmiiller, J. 
Chem. Phys. 116, 897 (2002). 

[17] A. F. Voter, Phys. Rev. Lett. 78, 3908 (1997) 

[18] A. F. Voter, J. Chem. Phys. 106, 4665 (1997). 

[19] D. Hamelberg, J. Mongan, and J. A. McCammon, J. 
Chem. Phys. 120 , 11919 (2004). 

[20] P. Tiwary and A. van de Walle, Phys. Rev. B 84, 100301 

( 2011 ) 

[21] P. Tiwary and A. van de Walle, Phys. Rev. B 87, 094304 
(2013) 

[22] P. Tiwary and M. Parrinello, Phys. Rev. Lett. Ill, 
230602 (2013) 

[23] R. A. Miron and K. A. Fichthorn, J. Chem. Phys. 119, 
6210 (2003). 

[24] P. G. Bolhuis, D. Chandler, C. Dellago, and P. L. 
Geissler, Annu. Rev. Phys. Chem. 53, 291 (2002). 

[25] A. K. Faradjian and R. Fiber, J. Chem. Phys. 120, 10880 
(2004). 

[26] T.-Q. Yu, A. Bucci, E. Vanden-Eijnden, and C. Abrams, 
Biophys. J. 108, 216a (2015). 

[27] X.-M. Bai, A. F. Voter, R. G. Hoagland, M. Nastasi, and 
B. P. Uberuaga, Science 327, 1631 (2010). 

[28] L. Maragliano, G. Cottone, G. Ciccotti, and E. Vanden- 
Eijnden, J. Am. Chem. Soc. 132, 1010 (2009). 

[29] A. F. Voter, F. Montalenti, and T. C. Germann, Ann. 
Rev. Mater. Res. 32, 321 (2002). 

[30] G. Henkelman and H. Jonsson, J. Chem. Phys. 115, 9657 

( 2001 ). 

[31] M. Salvalaglio, P. Tiwary, and M. Parrinello, J. Chem. 
Theory Comput. 10 , 1420 (2014). 

[32] P. Tiwary, V. Limongelli, M. Salvalaglio, and M. Par¬ 
rinello, Proc. Natl. Acad. Sci. U.S.A. , 201424461 (2015). 

[33] M. U. Bohner, J. Zeman, J. Smiatek, A. Arnold, and 
J. Kastner, J. Chem. Phys. 140, 074109 (2014). 

[34] J. Schneider and K. Reuter, J. Phys. Chem. Lett. 5, 3859 
(2014). 

[35] F. Sicard, N. Destainville, and M. Manghi, J. Chem. 
Phys. 142, 034903 (2015). 

[36] O. Valsson and M. Parrinello, Phys. Rev. Lett. 113, 
090601 (2014). 

[37] M. M. Steiner, P.-A. Genilloud, and J. W. Wilkins, Phys. 
Rev. B 57, 10236 (1998). 


















6 


[38] A. Chaimovich and M. S. Shell, J. Chem. Phys. 134, 
094112 (2011). 

[39] 1. Bilionis and P.-S. Koutsourelakis, J. Comput. Phys. 
231, 3849 (2012). 

[40] O. Valsson and M. Parrinello, J. Chem. Theory Comput. 
11, 1996 (2015). 


[41] F. Bach and E. Moulines, in Advances in Neural Informa¬ 
tion Processing Systems 26, edited by C. Burges, L. Bot- 
tou, M. Welling, Z. Ghahramani, and K. Weinberger 
(Curran Associates, Inc., Red Hook, NY, 2013) pp. 773- 
781. 

[42] G. A. Tribello, M. Bonomi, D. Branduardi, C. Camilloni, 
and G. Bussi, Comput. Phys. Commun. 185, 604 (2014). 



